knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 10'000 per chain to achieve 40'000
warmup = 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,24)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,24+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,47)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,10+5),
  "Conditional Between-Person Effects" = c(11+5,17+5),
  
  "Hurdle Within-Person Effects" = c(18+5,25+5),
  "Hurdle Between-Person Effects" = c(26+5,32+5),
  
  "Random Effects" = c(33+5, 46+5), 
  "Additional Parameters" = c(47+5,47+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_NOAR", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 40000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -10490.4 166.4
## p_loo       181.3   6.2
## looic     20980.8 332.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 3736, p-value = 0.8
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005353319 0.0032119914
## sample estimates:
## outlier frequency (expected: 0.00166220556745182 ) 
##                                        0.001873662
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd', 'ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, rope_range = rope_range_continuous,
## model_rows_fixed = model_rows_fixed_hu, : Coefficients were exponentiated.
## Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 47.95*** 3.03 [42.23, 54.34] 100.00% [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.001 7325 14978
Hurdle Intercept 0.85 0.14 [ 0.61, 1.18] 83.80% [0.84, 1.20] 0.515 0.099 Strong Evidence for Null 1.000 4644 10721
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.03 [ 0.97, 1.08] 83.24% [0.92, 1.08] 0.969 0.017 Very Strong Evidence for Null 1.000 11179 19513
Daily persuasion utilized (partner’s view) 1.03 0.02 [ 0.98, 1.08] 90.16% [0.92, 1.08] 0.976 0.023 Very Strong Evidence for Null 1.000 16280 22705
Daily pressure experienced 0.89* 0.04 [ 0.80, 0.98] 98.71% [0.92, 1.08] 0.215 0.281 Moderate Evidence for Null 1.000 25549 25127
Daily pressure utilized (partner’s view) 0.94 0.04 [ 0.86, 1.03] 91.68% [0.92, 1.08] 0.650 0.045 Strong Evidence for Null 1.000 27538 22984
Daily pushing experienced 1.03 0.04 [ 0.96, 1.10] 76.93% [0.92, 1.08] 0.933 0.019 Very Strong Evidence for Null 1.000 22879 25906
Daily pushing utilized (partner’s view) 0.99 0.03 [ 0.93, 1.05] 63.65% [0.92, 1.08] 0.978 0.013 Very Strong Evidence for Null 1.000 22162 25742
Day 1.01 0.06 [ 0.89, 1.14] 54.95% [0.92, 1.08] 0.789 0.025 Very Strong Evidence for Null 1.000 39583 31325
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.01 0.16 [ 0.74, 1.39] 53.40% [0.92, 1.08] 0.384 0.063 Strong Evidence for Null 1.000 5679 11615
Mean persuasion utilized (partner’s view) 0.98 0.15 [ 0.72, 1.35] 53.95% [0.92, 1.08] 0.384 0.064 Strong Evidence for Null 1.001 5613 11737
Mean pressure experienced 1.14 0.21 [ 0.80, 1.64] 76.54% [0.92, 1.08] 0.263 0.093 Strong Evidence for Null 1.001 8538 16974
Mean pressure utilized (partner’s view) 0.88 0.17 [ 0.61, 1.29] 74.45% [0.92, 1.08] 0.261 0.094 Strong Evidence for Null 1.000 8892 17039
Mean pushing experienced 1.33 0.31 [ 0.84, 2.09] 88.55% [0.92, 1.08] 0.128 0.195 Moderate Evidence for Null 1.001 8468 16152
Mean pushing utilized (partner’s view) 1.40 0.33 [ 0.88, 2.24] 92.27% [0.92, 1.08] 0.094 0.265 Moderate Evidence for Null 1.001 7823 16081
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.53*** 0.10 [ 1.36, 1.75] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 18293 23098
Hu Daily persuasion utilized (partner’s view) 1.32*** 0.08 [ 1.19, 1.50] 100.00% [0.84, 1.20] 0.038 >100 Overwhelming Evidence 1.000 21760 23402
Hu Daily pressure experienced 0.82 0.13 [ 0.58, 1.13] 89.75% [0.84, 1.20] 0.431 0.183 Moderate Evidence for Null 1.000 23948 22959
Hu Daily pressure utilized (partner’s view) 1.47* 0.27 [ 1.05, 2.32] 98.77% [0.84, 1.20] 0.115 1.184 Weak Evidence 1.000 22679 16886
Hu Daily pushing experienced 1.71*** 0.28 [ 1.27, 2.46] 99.97% [0.84, 1.20] 0.009 25.585 Strong Evidence 1.000 18260 21188
Hu Daily pushing utilized (partner’s view) 1.83*** 0.23 [ 1.46, 2.44] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 21693 23686
Hu Day 0.92 0.12 [ 0.71, 1.19] 74.29% [0.84, 1.20] 0.739 0.081 Strong Evidence for Null 1.000 42026 30353
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.19 0.46 [ 0.54, 2.60] 67.33% [0.84, 1.20] 0.321 0.217 Moderate Evidence for Null 1.001 4572 9486
Hu Mean persuasion utilized (partner’s view) 1.18 0.46 [ 0.53, 2.56] 66.74% [0.84, 1.20] 0.327 0.213 Moderate Evidence for Null 1.001 4614 9826
Hu Mean pressure experienced 0.31** 0.13 [ 0.13, 0.74] 99.63% [0.84, 1.20] 0.011 7.828 Moderate Evidence 1.000 7034 15087
Hu Mean pressure utilized (partner’s view) 0.57 0.25 [ 0.24, 1.36] 90.16% [0.84, 1.20] 0.143 0.504 Weak Evidence for Null 1.000 7037 15000
Hu Mean pushing experienced 2.70 1.53 [ 0.87, 8.16] 95.77% [0.84, 1.20] 0.056 1.372 Weak Evidence 1.000 7198 13867
Hu Mean pushing utilized (partner’s view) 2.79 1.57 [ 0.89, 8.46] 96.16% [0.84, 1.20] 0.051 1.415 Weak Evidence 1.000 6943 12974
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.32 0.04 [0.24, 0.42] NA NA NA NA NA 1.000 8803 15643
sd(Hurdle Intercept) 0.89 0.12 [0.69, 1.18] NA NA NA NA NA 1.000 7673 14307
sd(Daily persuasion experienced) 0.12 0.02 [0.08, 0.17] NA NA NA NA NA 1.000 16113 24856
sd(Daily persuasion utilized (partner’s view)) 0.09 0.02 [0.05, 0.13] NA NA NA NA NA 1.000 19267 21079
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.001 12569 14433
sd(Daily pressure utilized (partner’s view)) 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 15590 14458
sd(Daily pushing experienced) 0.11 0.04 [0.04, 0.19] NA NA NA NA NA 1.000 11320 7606
sd(Daily pushing utilized (partner’s view)) 0.09 0.03 [0.02, 0.17] NA NA NA NA NA 1.000 10989 6835
sd(Hu Daily persuasion experienced) 0.18 0.08 [0.02, 0.34] NA NA NA NA NA 1.000 8938 6082
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 0.08 [0.02, 0.33] NA NA NA NA NA 1.000 9475 7735
sd(Hu Daily pressure experienced) 0.25 0.21 [0.01, 0.85] NA NA NA NA NA 1.001 11940 15619
sd(Hu Daily pressure utilized (partner’s view)) 0.28 0.25 [0.01, 1.00] NA NA NA NA NA 1.000 10937 15282
sd(Hu Daily pushing experienced) 0.62 0.18 [0.32, 1.08] NA NA NA NA NA 1.000 15379 20817
sd(Hu Daily pushing utilized (partner’s view)) 0.31 0.15 [0.04, 0.64] NA NA NA NA NA 1.000 10681 7340
Additional Parameters
sigma 0.68 0.01 [0.66, 0.71] NA NA NA NA NA 1.000 38393 29400 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.83). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 158 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 158 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00654

$persuasion_partner_cw

## Warning: Removed 137 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 129 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00537

$pressure_self_cw

## Warning: Removed 45 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 51 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00939

$pressure_partner_cw

## Warning: Removed 45 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 87 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0156

$pushing_self_cw

## Warning: Removed 35 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 145 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0153

$pushing_partner_cw

## Warning: Removed 26 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 63 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0125

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  rmarkdown::render(
    "C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/BeautifulPlotWithNote.Rmd", 
    output_file = paste0('C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/Graphic_', effname, '.pdf'),
    params = list(
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}
marginaleffects::avg_predictions(pa_sub)
## 
##  Estimate 2.5 % 97.5 %
##      30.9  29.4   32.5
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 40000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -18781.8  68.4
## p_loo        92.1   4.5
## looic     37563.7 136.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 24, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001341025 0.004652382
## sample estimates:
## outlier frequency (expected: 0.00287383877734492 ) 
##                                        0.007192089
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, rope_range = rope_range_log,
## model_rows_fixed = model_rows_fixed, : Coefficients were exponentiated. Double
## check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 117.27*** 6.24 [105.41, 130.66] 100.00% [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 5582 10902
Within-Person Effects
Daily persuasion experienced 1.03 0.02 [ 1.00, 1.06] 96.53% [0.94, 1.07] 0.990 0.033 Strong Evidence for Null 1.000 22787 26727
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.05] 88.83% [0.94, 1.07] 0.997 0.014 Very Strong Evidence for Null 1.000 28722 29910
Daily pressure experienced 0.94 0.03 [ 0.88, 1.01] 95.99% [0.94, 1.07] 0.547 0.066 Strong Evidence for Null 1.000 41819 28808
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.92, 1.05] 71.66% [0.94, 1.07] 0.905 0.016 Very Strong Evidence for Null 1.000 46892 31087
Daily pushing experienced 1.03 0.03 [ 0.98, 1.08] 90.00% [0.94, 1.07] 0.907 0.025 Very Strong Evidence for Null 1.000 32795 27546
Daily pushing utilized (partner’s view) 1.02 0.02 [ 0.97, 1.06] 77.40% [0.94, 1.07] 0.986 0.011 Very Strong Evidence for Null 1.000 46123 31056
Day 0.97 0.03 [ 0.91, 1.04] 78.37% [0.94, 1.07] 0.849 0.019 Very Strong Evidence for Null 1.000 63974 30299
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 100.00% [0.94, 1.07] 1.000 0.106 Moderate Evidence for Null 1.000 41200 26262
Between-Person Effects
Mean persuasion experienced 1.10 0.16 [ 0.83, 1.47] 75.59% [0.94, 1.07] 0.275 0.074 Strong Evidence for Null 1.001 5262 10215
Mean persuasion utilized (partner’s view) 0.98 0.14 [ 0.73, 1.31] 55.00% [0.94, 1.07] 0.349 0.057 Strong Evidence for Null 1.001 5260 10280
Mean pressure experienced 0.98 0.14 [ 0.73, 1.31] 56.39% [0.94, 1.07] 0.339 0.059 Strong Evidence for Null 1.001 7201 14024
Mean pressure utilized (partner’s view) 0.96 0.14 [ 0.72, 1.28] 60.24% [0.94, 1.07] 0.339 0.059 Strong Evidence for Null 1.001 6726 13627
Mean pushing experienced 0.97 0.20 [ 0.64, 1.47] 56.15% [0.94, 1.07] 0.242 0.083 Strong Evidence for Null 1.000 7856 13134
Mean pushing utilized (partner’s view) 1.24 0.25 [ 0.83, 1.88] 85.86% [0.94, 1.07] 0.141 0.147 Moderate Evidence for Null 1.000 7770 12879
Mean weartime 1.00 0.00 [ 1.00, 1.00] 91.16% [0.94, 1.07] 1.000 0.000 Very Strong Evidence for Null 1.000 50538 35353
Random Effects
sd(Intercept) 0.30 0.04 [0.24, 0.40] NA NA NA NA NA 1.000 7902 14886
sd(Daily persuasion experienced) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 20801 17210
sd(Daily persuasion utilized (partner’s view)) 0.06 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 20659 19450
sd(Daily pressure experienced) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 17044 16600
sd(Daily pressure utilized (partner’s view)) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 23038 19677
sd(Daily pushing experienced) 0.08 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 8952 10098
sd(Daily pushing utilized (partner’s view)) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 13943 18163
Additional Parameters
sigma 0.57 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 60296 29184 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.89), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.72), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.78), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.83). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5186.2  59.2
## p_loo        75.4   3.3
## looic     10372.5 118.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  5780    
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 27, observations = 3736, p-value =
## 2.533e-08
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.004767871 0.010497581
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.007226981
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summarize_brms(
  mood_gauss, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.70*** 0.10 [ 3.50, 3.91] 100.00% [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.002 4495 9545
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.04, 0.05] 54.54% [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 35462 27914
Daily persuasion utilized (partner’s view) 0.02 0.02 [-0.02, 0.07] 82.80% [-0.11, 0.11] 1.000 0.007 Very Strong Evidence for Null 1.000 29416 30178
Daily pressure experienced -0.04 0.05 [-0.15, 0.07] 76.90% [-0.11, 0.11] 0.919 0.013 Very Strong Evidence for Null 1.000 35865 27332
Daily pressure utilized (partner’s view) -0.02 0.05 [-0.14, 0.08] 66.80% [-0.11, 0.11] 0.934 0.012 Very Strong Evidence for Null 1.000 33073 24399
Daily pushing experienced 0.02 0.03 [-0.04, 0.08] 75.78% [-0.11, 0.11] 0.997 0.008 Very Strong Evidence for Null 1.000 41974 29904
Daily pushing utilized (partner’s view) 0.08* 0.03 [ 0.01, 0.14] 98.45% [-0.11, 0.11] 0.867 0.080 Strong Evidence for Null 1.000 31271 28932
Day 0.26*** 0.06 [ 0.15, 0.37] 100.00% [-0.11, 0.11] 0.005 >100 Overwhelming Evidence 1.000 59723 29576
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.35 0.28 [-0.20, 0.89] 89.30% [-0.11, 0.11] 0.152 0.123 Moderate Evidence for Null 1.000 4262 9105
Mean persuasion utilized (partner’s view) 0.23 0.28 [-0.31, 0.78] 79.91% [-0.11, 0.11] 0.227 0.080 Strong Evidence for Null 1.000 4267 8956
Mean pressure experienced -0.32 0.27 [-0.86, 0.23] 87.62% [-0.11, 0.11] 0.170 0.108 Moderate Evidence for Null 1.000 4747 11106
Mean pressure utilized (partner’s view) -0.32 0.28 [-0.86, 0.23] 87.67% [-0.11, 0.11] 0.173 0.110 Moderate Evidence for Null 1.000 4951 11430
Mean pushing experienced 0.22 0.38 [-0.55, 1.01] 71.39% [-0.11, 0.11] 0.200 0.089 Strong Evidence for Null 1.001 6400 11111
Mean pushing utilized (partner’s view) 0.35 0.38 [-0.41, 1.14] 82.26% [-0.11, 0.11] 0.156 0.118 Moderate Evidence for Null 1.001 6362 11744
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.08 [0.48, 0.78] NA NA NA NA NA 1.001 7483 13998
sd(Daily persuasion experienced) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.001 10758 15861
sd(Daily persuasion utilized (partner’s view)) 0.08 0.03 [0.02, 0.13] NA NA NA NA NA 1.000 11761 9894
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 16544 18107
sd(Daily pressure utilized (partner’s view)) 0.08 0.07 [0.00, 0.26] NA NA NA NA NA 1.000 14798 17566
sd(Daily pushing experienced) 0.05 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 13961 15345
sd(Daily pushing utilized (partner’s view)) 0.07 0.05 [0.00, 0.17] NA NA NA NA NA 1.000 13129 14428
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 63319 29501 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.82), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.82), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.89). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw 

marginaleffects::avg_predictions(mood_gauss)
## 
##  Estimate 2.5 % 97.5 %
##      3.82  3.79   3.85
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_NOARNOAR_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -682.1 31.9
## p_loo        73.7  5.4
## looic      1364.3 63.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     751   99.3%   903     
##    (0.7, 1]   (bad)        5    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.04
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.001322751
## sample estimates:
## outlier frequency (expected: 0.000436507936507937 ) 
##                                         0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summarize_brms(
  reactance_ordinal, 
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.84*** 0.98 [ 2.34, 6.49] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 28792 29901
Intercept[2] 8.32*** 2.24 [ 4.96, 14.45] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 29900 30020
Intercept[3] 23.17*** 6.80 [ 13.17, 42.42] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 31532 31530
Intercept[4] 101.19*** 35.37 [ 52.40, 208.53] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 34784 31711
Intercept[5] 3422.14*** 2190.68 [1091.76, 13601.17] 100.00% [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 47511 31296
Within-Person Effects
Daily persuasion experienced 0.85* 0.07 [ 0.71, 0.99] 97.97% [0.84, 1.20] 0.568 0.277 Moderate Evidence for Null 1.000 33454 27707
Daily persuasion utilized (partner’s view) 1.03 0.10 [ 0.84, 1.24] 60.38% [0.84, 1.20] 0.922 0.040 Strong Evidence for Null 1.000 30268 27064
Daily pressure experienced 1.86* 0.36 [ 1.18, 2.69] 99.38% [0.84, 1.20] 0.027 3.012 Moderate Evidence 1.000 19636 23573
Daily pressure utilized (partner’s view) 1.23 0.29 [ 0.70, 2.08] 80.82% [0.84, 1.20] 0.385 0.149 Moderate Evidence for Null 1.000 25256 21310
Daily pushing experienced 1.16 0.11 [ 0.96, 1.42] 94.52% [0.84, 1.20] 0.618 0.136 Moderate Evidence for Null 1.000 28471 27169
Daily pushing utilized (partner’s view) 0.92 0.11 [ 0.71, 1.17] 76.80% [0.84, 1.20] 0.758 0.062 Strong Evidence for Null 1.000 32368 25496
Day 1.47 0.51 [ 0.77, 2.86] 87.16% [0.84, 1.20] 0.227 0.260 Moderate Evidence for Null 1.000 45161 31117
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.12 0.57 [ 0.41, 3.14] 59.12% [0.84, 1.20] 0.268 0.212 Moderate Evidence for Null 1.000 13321 21398
Mean persuasion utilized (partner’s view) 1.38 0.79 [ 0.45, 4.35] 71.27% [0.84, 1.20] 0.216 0.263 Moderate Evidence for Null 1.000 14141 21745
Mean pressure experienced 3.48* 1.90 [ 1.18, 10.66] 98.83% [0.84, 1.20] 0.022 3.043 Moderate Evidence 1.000 15342 23058
Mean pressure utilized (partner’s view) 1.18 0.67 [ 0.37, 3.58] 61.23% [0.84, 1.20] 0.236 0.241 Moderate Evidence for Null 1.000 15692 22489
Mean pushing experienced 1.20 0.89 [ 0.28, 5.59] 59.92% [0.84, 1.20] 0.189 0.305 Weak Evidence for Null 1.000 18260 25276
Mean pushing utilized (partner’s view) 0.11* 0.10 [ 0.02, 0.64] 99.30% [0.84, 1.20] 0.008 8.595 Moderate Evidence 1.000 23596 28667
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.81 0.20 [0.48, 1.26] NA NA NA NA NA 1.000 14049 23078
sd(Daily persuasion experienced) 0.16 0.12 [0.01, 0.42] NA NA NA NA NA 1.000 7067 14443
sd(Daily persuasion utilized (partner’s view)) 0.21 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 10643 14404
sd(Daily pressure experienced) 0.55 0.25 [0.10, 1.14] NA NA NA NA NA 1.001 9501 8812
sd(Daily pressure utilized (partner’s view)) 0.40 0.38 [0.02, 1.55] NA NA NA NA NA 1.000 10032 16731
sd(Daily pushing experienced) 0.20 0.13 [0.01, 0.51] NA NA NA NA NA 1.000 10833 13309
sd(Daily pushing utilized (partner’s view)) 0.14 0.13 [0.01, 0.59] NA NA NA NA NA 1.000 15145 18957
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.76), b_Intercept[4] and b_Intercept[3] (r = 0.83),
##   b_pressure_self_cb and b_persuasion_self_cb (r = 0.72),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.79). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModels_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw 

marginaleffects::avg_predictions(reactance_ordinal)
## 
##  Group Estimate   2.5 % 97.5 %
##      0  0.68193 0.65457 0.7076
##      1  0.09412 0.07508 0.1157
##      2  0.08361 0.06612 0.1041
##      3  0.06984 0.05464 0.0874
##      4  0.06395 0.05130 0.0775
##      5  0.00527 0.00171 0.0115
## 
## Type:  response 
## Columns: group, estimate, conf.low, conf.high

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_NOARNOAR_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
check_brms(is_reactance)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 34 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -363.6 16.0
## p_loo        80.3  6.0
## looic       727.2 32.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     722   95.5%   476     
##    (0.7, 1]   (bad)       33    4.4%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.001322751
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.29*** 0.08 [0.16, 0.50] 100.00% [0.83, 1.20] 0.000 >100 Overwhelming Evidence 1.000 43132 34475
Within-Person Effects
Daily persuasion experienced 0.84 0.08 [0.69, 1.01] 96.47% [0.83, 1.20] 0.552 0.191 Moderate Evidence for Null 1.000 42037 30341
Daily persuasion utilized (partner’s view) 1.12 0.16 [0.85, 1.53] 78.62% [0.83, 1.20] 0.665 0.076 Strong Evidence for Null 1.000 34121 28236
Daily pressure experienced 1.99* 0.63 [1.03, 4.47] 97.91% [0.83, 1.20] 0.050 1.446 Weak Evidence 1.000 28128 23067
Daily pressure utilized (partner’s view) 1.41 0.59 [0.59, 4.18] 80.24% [0.83, 1.20] 0.241 0.239 Moderate Evidence for Null 1.000 27770 24147
Daily pushing experienced 1.27* 0.15 [1.01, 1.63] 97.89% [0.83, 1.20] 0.309 0.391 Weak Evidence for Null 1.000 45536 30516
Daily pushing utilized (partner’s view) 0.89 0.17 [0.60, 1.31] 73.50% [0.83, 1.20] 0.570 0.093 Strong Evidence for Null 1.000 43340 30983
Day 1.64 0.63 [0.76, 3.51] 89.82% [0.83, 1.20] 0.165 0.350 Weak Evidence for Null 1.000 66504 30884
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.99 1.20 [0.60, 6.84] 87.26% [0.83, 1.20] 0.127 0.462 Weak Evidence for Null 1.000 22054 28793
Mean persuasion utilized (partner’s view) 1.88 1.25 [0.52, 7.23] 82.88% [0.83, 1.20] 0.137 0.419 Weak Evidence for Null 1.000 23180 27564
Mean pressure experienced 17.42** 18.55 [2.39, 155.22] 99.77% [0.83, 1.20] 0.003 22.823 Strong Evidence 1.000 30159 28196
Mean pressure utilized (partner’s view) 2.32 2.55 [0.25, 19.47] 77.24% [0.83, 1.20] 0.099 0.593 Weak Evidence for Null 1.000 26009 28694
Mean pushing experienced 0.82 0.81 [0.12, 6.10] 58.04% [0.83, 1.20] 0.144 0.408 Weak Evidence for Null 1.000 28811 29362
Mean pushing utilized (partner’s view) 0.08* 0.08 [0.01, 0.65] 99.01% [0.83, 1.20] 0.009 7.113 Moderate Evidence 1.000 32814 31176
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.15 0.25 [0.74, 1.74] NA NA NA NA NA 1.000 17091 25690
sd(Daily persuasion experienced) 0.20 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 9903 17543
sd(Daily persuasion utilized (partner’s view)) 0.48 0.20 [0.12, 0.97] NA NA NA NA NA 1.000 14070 15231
sd(Daily pressure experienced) 1.06 0.55 [0.14, 2.41] NA NA NA NA NA 1.000 9656 11046
sd(Daily pressure utilized (partner’s view)) 0.81 0.67 [0.04, 2.73] NA NA NA NA NA 1.000 14835 20412
sd(Daily pushing experienced) 0.23 0.16 [0.01, 0.59] NA NA NA NA NA 1.000 14895 17689
sd(Daily pushing utilized (partner’s view)) 0.25 0.22 [0.01, 0.93] NA NA NA NA NA 1.001 17388 22176
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModels_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

marginaleffects::avg_predictions(is_reactance)
## 
##  Estimate 2.5 % 97.5 %
##     0.329 0.304  0.355
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.47      0.39    -0.15     1.13       8.81
##   Post.Prob Star
## 1       0.9     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsFinal.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 47.95*** [42.23, 54.34] 100.00% 117.27*** [105.41, 130.66] 100.00% 3.70*** [ 3.50, 3.91] 100.00% NA NA NA 0.29*** [0.16, 0.50] 100.00%
Hurdle Intercept 0.85 [ 0.61, 1.18] 83.80% NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.03 [ 0.97, 1.08] 83.24% 1.03 [ 1.00, 1.06] 96.53% 0.00 [-0.04, 0.05] 54.54% 0.85* [ 0.71, 0.99] 97.97% 0.84 [0.69, 1.01] 96.47%
Daily persuasion utilized (partner’s view) 1.03 [ 0.98, 1.08] 90.16% 1.02 [ 0.99, 1.05] 88.83% 0.02 [-0.02, 0.07] 82.80% 1.03 [ 0.84, 1.24] 60.38% 1.12 [0.85, 1.53] 78.62%
Daily pressure experienced 0.89* [ 0.80, 0.98] 98.71% 0.94 [ 0.88, 1.01] 95.99% -0.04 [-0.15, 0.07] 76.90% 1.86* [ 1.18, 2.69] 99.38% 1.99* [1.03, 4.47] 97.91%
Daily pressure utilized (partner’s view) 0.94 [ 0.86, 1.03] 91.68% 0.98 [ 0.92, 1.05] 71.66% -0.02 [-0.14, 0.08] 66.80% 1.23 [ 0.70, 2.08] 80.82% 1.41 [0.59, 4.18] 80.24%
Daily pushing experienced 1.03 [ 0.96, 1.10] 76.93% 1.03 [ 0.98, 1.08] 90.00% 0.02 [-0.04, 0.08] 75.78% 1.16 [ 0.96, 1.42] 94.52% 1.27* [1.01, 1.63] 97.89%
Daily pushing utilized (partner’s view) 0.99 [ 0.93, 1.05] 63.65% 1.02 [ 0.97, 1.06] 77.40% 0.08* [ 0.01, 0.14] 98.45% 0.92 [ 0.71, 1.17] 76.80% 0.89 [0.60, 1.31] 73.50%
Day 1.01 [ 0.89, 1.14] 54.95% 0.97 [ 0.91, 1.04] 78.37% 0.26*** [ 0.15, 0.37] 100.00% 1.47 [ 0.77, 2.86] 87.16% 1.64 [0.76, 3.51] 89.82%
Daily weartime NA NA NA 1.00*** [ 1.00, 1.00] 100.00% NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.01 [ 0.74, 1.39] 53.40% 1.10 [ 0.83, 1.47] 75.59% 0.35 [-0.20, 0.89] 89.30% 1.12 [ 0.41, 3.14] 59.12% 1.99 [0.60, 6.84] 87.26%
Mean persuasion utilized (partner’s view) 0.98 [ 0.72, 1.35] 53.95% 0.98 [ 0.73, 1.31] 55.00% 0.23 [-0.31, 0.78] 79.91% 1.38 [ 0.45, 4.35] 71.27% 1.88 [0.52, 7.23] 82.88%
Mean pressure experienced 1.14 [ 0.80, 1.64] 76.54% 0.98 [ 0.73, 1.31] 56.39% -0.32 [-0.86, 0.23] 87.62% 3.48* [ 1.18, 10.66] 98.83% 17.42** [2.39, 155.22] 99.77%
Mean pressure utilized (partner’s view) 0.88 [ 0.61, 1.29] 74.45% 0.96 [ 0.72, 1.28] 60.24% -0.32 [-0.86, 0.23] 87.67% 1.18 [ 0.37, 3.58] 61.23% 2.32 [0.25, 19.47] 77.24%
Mean pushing experienced 1.33 [ 0.84, 2.09] 88.55% 0.97 [ 0.64, 1.47] 56.15% 0.22 [-0.55, 1.01] 71.39% 1.20 [ 0.28, 5.59] 59.92% 0.82 [0.12, 6.10] 58.04%
Mean pushing utilized (partner’s view) 1.40 [ 0.88, 2.24] 92.27% 1.24 [ 0.83, 1.88] 85.86% 0.35 [-0.41, 1.14] 82.26% 0.11* [ 0.02, 0.64] 99.30% 0.08* [0.01, 0.65] 99.01%
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 91.16% NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.53*** [ 1.36, 1.75] 100.00% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.32*** [ 1.19, 1.50] 100.00% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 0.82 [ 0.58, 1.13] 89.75% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.47* [ 1.05, 2.32] 98.77% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 1.71*** [ 1.27, 2.46] 99.97% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 1.83*** [ 1.46, 2.44] 100.00% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.92 [ 0.71, 1.19] 74.29% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.19 [ 0.54, 2.60] 67.33% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.18 [ 0.53, 2.56] 66.74% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.31** [ 0.13, 0.74] 99.63% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.57 [ 0.24, 1.36] 90.16% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 2.70 [ 0.87, 8.16] 95.77% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 2.79 [ 0.89, 8.46] 96.16% NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.32 [0.24, 0.42] NA 0.30 [0.24, 0.40] NA 0.60 [0.48, 0.78] NA 0.81 [0.48, 1.26] NA 1.15 [0.74, 1.74] NA
sd(Hurdle Intercept) 0.89 [0.69, 1.18] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.12 [0.08, 0.17] NA 0.05 [0.02, 0.08] NA 0.04 [0.00, 0.10] NA 0.16 [0.01, 0.42] NA 0.20 [0.01, 0.51] NA
sd(Daily persuasion utilized (partner’s view)) 0.09 [0.05, 0.13] NA 0.06 [0.03, 0.09] NA 0.08 [0.02, 0.13] NA 0.21 [0.01, 0.51] NA 0.48 [0.12, 0.97] NA
sd(Daily pressure experienced) 0.07 [0.00, 0.24] NA 0.04 [0.00, 0.13] NA 0.07 [0.00, 0.24] NA 0.55 [0.10, 1.14] NA 1.06 [0.14, 2.41] NA
sd(Daily pressure utilized (partner’s view)) 0.06 [0.00, 0.18] NA 0.04 [0.00, 0.11] NA 0.08 [0.00, 0.26] NA 0.40 [0.02, 1.55] NA 0.81 [0.04, 2.73] NA
sd(Daily pushing experienced) 0.11 [0.04, 0.19] NA 0.08 [0.01, 0.15] NA 0.05 [0.00, 0.14] NA 0.20 [0.01, 0.51] NA 0.23 [0.01, 0.59] NA
sd(Daily pushing utilized (partner’s view)) 0.09 [0.02, 0.17] NA 0.04 [0.00, 0.10] NA 0.07 [0.00, 0.17] NA 0.14 [0.01, 0.59] NA 0.25 [0.01, 0.93] NA
sd(Hu Daily persuasion experienced) 0.18 [0.02, 0.34] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 [0.02, 0.33] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.25 [0.01, 0.85] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.28 [0.01, 1.00] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.62 [0.32, 1.08] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.31 [0.04, 0.64] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.68 [0.66, 0.71] NA 0.57 [0.56, 0.59] NA 0.96 [0.94, 0.98] NA NA NA NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()